# rm(list = ls())

gc()

# source("00_pckgs.R")

dir.create("data/mc", recursive = TRUE)

# set.seed(1917)

for (j in 1:1000){
  
  # sample correlation in ideal points
  rho <- sample(seq(-.8, .8, by = .1), 1)
  
  # sample weight of dimension
  taus <- sample(seq(.15, .45, by = .1), 1)
  tau <- sort(c(1 - taus, taus), T)
  
  # sample level of nonsep
  dl <- sample(seq(-.8, .8, by = .1), 1)
  
  # covariates for ideal points
  id_i <- MASS::mvrnorm(100, c(0, 0), matrix(c(1, rho, rho, 1), ncol = 2))
  id_i <- apply(id_i, 2, function(x){
    ifelse(x < quantile(x, 1/3), -1,
           ifelse(x > quantile(x, 2/3), 1,
                  0))
  })
  
  # ideal points based on covariates
  # roughly 2/3 systematic, 1/3 stochastic
  id <- cbind(rnorm(100, id_i[,1], sqrt(1/3)),
              rnorm(100, id_i[,2], sqrt(1/3)))
  
  # disc and diff pars of roll call votes
  prop <- MASS::mvrnorm(250, c(0, 0), matrix(c(1, 0, 0, 1), ncol = 2))
  prop_diff <- rnorm(250) * 2
  
  # standardize
  id <- apply(id, 2, function(x) (x - mean(x)) / sd(x))
  prop <- apply(prop, 2, function(x) (x - mean(x)) / sd(x))
  
  # scale weight matrix so that systematic component accounts for roughly 2/3 of variance in resulting ystar; var(ystar) = var(y) + var(error)
  w <- diag(sqrt(tau)) %*% matrix(c(1, dl, dl, 1), ncol = 2) %*% diag(sqrt(tau))
  w_sc <- w / sqrt(sum(w**2)) * sqrt(2)
  
  y <- id %*% w_sc %*% t(prop)
  
  for (i in 1:250){
    y[,i] <- y[,i] + prop_diff[i] + rnorm(100)
  }
  
  y_b <- pnorm(y)
  
  for (k in 1:100){
    for (i in 1:250){
      y_b[k,i] <- rbinom(1, 1, y_b[k,i])
    }
  }
  save(dl, taus, rho,
       prop, prop_diff, id, id_i,
       y_b,
       file = paste0("data/mc/", j, "_run.Rda"))
}
